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We have studied the soliton propagation through a segment containing random point-like scatter- 
ers. In the limit of small concentration of scatterers when the mean distance between the scatterers 
is larger than the soliton width a method has been developed for obtaining the statistical charac- 
teristics of the soliton transmission through the segment. The method is applicable for any classical 
particle traversing through a disordered segment with the given velocity transformation after each 
act of scattering. In the case of weak scattering and relatively short disordered segment the trans- 
mission time delay of a fast soliton is mostly determined by the shifts of the soliton center after each 
act of scattering. For sufficiently long segments the main contribution to the delay is due to the 
shifts of the amplitude and velocity of a fast soliton after each scatterer. Corresponding crossover 
lengths for both cases of light and heavy solitons have been obtained. We have also calculated 
the exact probability density function of the soliton transmission time delay for a sufficiently long 
segment. In the case of weak identical scatterers the latter is a universal function which depends on 
a sole parameter - mean number of scatterers in a segment. 
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I. INTRODUCTION 



The interplay between disorder and nonlinearity has attracted the attention of physicists for more than twenty 
years [l|, 0] . The latest advancements in experimental techniques and computational facilities have brought about a 
new peak of interest in the problem and already led to a wider understanding of the corresponding phenomena [1| . 
Physical properties of disordered nonlinear systems reveal much more sensitivity to the details of the formulation of 
the problem than their linear counterparts. For example, in nonlinear systems there exist three different transmission 
regimes: beside the exponential decrease of the transmission coefficient with the length of the disordered segment 
H, Q (which is the sole possibility in the linear case) , power law decrease can occur [|| , and, in the case of sufficiently 
strong nonlinearity, the transmission coefficient may not even decrease at all [H, Q ■ 

One of the main concepts in the theory of nonlinear waves is a concept of a soliton 0-E| - a particle-like stable 
nonlinear excitation observed in many nonlinear physical systems. A particular example revealing the reciprocal 
action between the nonlinearity and disorder is a soliton propagation in random media. One can observe here an 
entire plethora of problems which differ from each other by (i) the ty pe of nonlinearity (e.g. nonlinear Schrodinger 
equation 

HEM!, sine-Gordon equation [IM1, so-called ^-systemMMIlJ), (ii) the way in which randomness 
enters the system ( e.g . in the random potential form (3, [ilj-flS . |2C|, dispersive terms [ll| , nonlinear terms 
[HI], external for ce fl9l0 i (iii) random structure of corresponding coefficients (e.g. randomly distributed localized 
impurities [i. [Tol [TH Il8l. |20| or finite-range inhomogeneities fljj - ll5l fl9l l2ll]). and (iv) the statistics of the disorder 
potential like e.g. white or colored Gaussian potential [l6[ or random step- wise process [12]. A noticeable share of 
the above mentioned papers is devoted to the soliton transmission through a ID disordered segment in the framework 
of nonlinear Schrodinger equation (NLSE). The latter is a ubiquitous nonlinear model appearing in many areas 
of contemporary physics. In the condensed matter physics NLSE occurs for example in the description of weakly 
nonlinear magnetization dynamics in ferromagnets with the easy-axis type anisotropy [12] ■ It is also one of the main 
models in the nonlinear fiber optics and the nonlinear fiber arrays [23l425| . A prominent features of the NLSE model 
is that in homogeneous (ordered) systems, this equation is completely integrable and possesses stable robust soliton 
solutions @, H[. 

Among the most common and widely explored sources of disorder are randomly distributed point-like scatterers 
which emulate the random short-scale imperfection of a media [2, [12, E3 ■ When the concentration of scatterers is 
small, that is to say the mean distance between scatterers is much larger than the size of a soliton (i.e. the case of 
"sparse scatterers" ) , the transmission of the soliton through a disordered segment can be considered as a sequence of 
individual events of passing through a single scattcrer. After each act of scattering, the transmitted soliton acquires 
an abrupt shift of its position as well as a change of both its energy and velocity (i.e. the scattering is generally 
inelastic). These changes influence the transmission characteristics of the soliton traversing through the total length 
of a sample, e.g. the time needed for the soliton to pass the entire segment (we call it a "transmission time" here). 
Assuming that the scattering is weak one can postulate that a soliton passing through a single scatterer, experiences 
only a slight change of its parameters, the measure of this "slightness" being the intensity of a scatterer. A single act 
of such a scattering was treated in detail in within the framework of perturbation theory [2] (see also Ref. (l3j 
where it was considered within another approach). When the perturbation method is applied, the resulting mapping 
mechanism for the soliton parameters is often called an "equivalent particle approach" due to the similarity of the 
resulting effect with the classical particle scattering [l4[ . 

In this paper we study the statistics of the soliton transmission time T through a ID disordered segment in the 
framework of NLSE. More precisely, we address the delay for the transmission time of an incident soliton, AT, 
occurring due to the presence of sparse random scatterers. Note, that often in the context of nonlinear optics the 
role of time is played by the spatial coordinate along the propagation of the optical beam, so for the so-called spatial 
solitons (23| this delay, AT, should be endowed with a different physical meaning. 

As mentioned above, the scattering itself is considered to be weak and a perturbative approach can be applied 
for the description of the soliton dynamics. In the first order with respect to weak scattering intensity, s, the time 
delay is related only to the shift in the soliton position after each scatterer. The corresponding total delay, AT, is 
then proportional to the number of scatterers within the segment (and hence the segment length L). On the other 
hand, the second order effects, i.e. those of the order of e 2 , bring about the change of the soliton velocity, and their 
corresponding contribution is proportional to L 2 . Therefore in the case of sufficiently long segment (the relevant 
inequalities will be given in the text), the latter contribution dominates, and the overall transmission delay should be 
calculated using the second order approximation with respect to the scatterer intensity. 

We develop here a simple but yet robust mathematical formalism allowing one to obtain the mean value, (AT), and 
variance, a 2 = ((AT) 2 ) — ((AT)) 2 , of the fluctuating transmission time delay. Some preliminary results concerning 
these first two moment have been published in [261 ] . But in the current paper not only we provide a detailed description 
of the method for the moments but also present the exact probability density, p(AT), and the probability P(AT < 
ATo) that the delay does not exceed a given value ATo. Such an approach has a certain merit on its own and can be 
used in the variety of similar random transmission problems where the velocity mapping between the two adjacent 
segments is given. It is pertinent to stress that in the current paper we consider the case of weak intensities of local 
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inhomogcncitics satisfying the inequality: e 2 n ~ e 2 L <C 1, where n is the typical number of scatterers. Thus the 
parameters of the transmitting soliton (i.e. its amplitude and velocity) change weakly, the measure of this weakness 
being the aforementioned parameter e 2 n. This case is somewhat opposite to that considered in Refs.fTH [Hj], where 
the limit e 2 L — > oo was implied so that the limiting values of the soliton parameters differ significantly from their 
initial values. 

The structure of the paper is as follows. We start with brief reminder of the basic properties of a NLSE soliton, 
consider soliton scattering on a single defect and recall the known expressions for the energy and the number of 
emitted quasiparticles (subsection III A[) used to obtain the formulas for velocity transformation (subsection IIIBI) . 
Our statistical model of randomly placed scatterers is introduced in section IIII1 where we describe the exact (non- 
pcrturbative) method for the calculation of various ensemble averages. This method is then applied in Section [IVI to 
the problem of soliton transmission through a disordered segment: we analyze dynamics of the soliton passing trough 
the segment (subsection IIV Ap , and obtain general expressions for the mean transmission time delay and its variance 
(subsection HVBp . The limiting case of weak scatterers is the subject of subsection IIV CI Here we obtain simple 
explicit formulas for the mean value and variance of the transmission time delay for both light and heavy solitons and 
estimate the length of the segment where second order contribution dominates over the first order one. The statistical 
properties of the delay time are studied in Section [Vl In its first part (subsection IV Ap we provide general formulas 
for various averages. The main result of this section is an explicit formula for the probability density function of 
transmission time delay obtained in subsection IV Bl In the last subsection IV CI we present the probability that delay 
does not exceed a given value. The results obtained are summarized in Conclusion. Appendices lAllDl contain the 
technical details of the calculations. 



II. SOLITON TRANSMISSION THROUGH A SINGLE WEAK SCATTERER 



A. General remarks. 



In this subsection we briefly recall some known results concerning the NLSE soliton transmission through the single 
weak point-like scatterer. Corresponding perturbed NLSE in the normalized dimensionless units reads [10(: 

iut + u xx + 2\u\ 2 u = ueS(x). — oo < x < oo. (1) 



Here u(x,t) is the complex field variable and the subscripts denote the partial derivatives with respect to time t and 
spatial coordinate x. The r.h.s. of this equation describes the influence of a single point scatterer with the intensity 
e placed at origin. Let us mention again that in some systems (like e.g. for the spatial solitons in planar waveguides) 
the physical meaning of independent variables t and x may be different. Also the meaning of the field variable u 
depends on the nature of the problem in hand: it can denote e.g. the deviation of the magnetization from the "easy" 
axis in ferromagnets or the envelope of the electromagnetic wave in optical fibers etc. 

The unperturbed NLSE, i.e. when we set e = in Eq.([T]), is completely integrable and possesses an infinite set 
of integrals of motion Q, including total energy E and the total number of bound waves (quasiparticles) N . The 
simplest one-soliton solution of NLSE reads as Q : 



cxp 



u s (x, t) = ia- 



vx 

Y 



t + ifio 



cosh [a(x + vt) — ipi] 



(2) 



and depends on four real parameters: a, v, ipo, (p±. We choose the soliton amplitude a (which also characterizes its 
spatial size), and velocity v to be positive. Such a selection corresponds to a soliton propagating in the negative 
direction of the x axis. 

Assuming the solution in the form of a single-soliton anstz, Eq. ([2]), the number of quasiparticles N s bounded in 
the soliton © and soliton energy E s are given by 



E„ = 



2a, 

v 2 2a 3 



(3) 
(4) 



In what follows we distinguish between the cases of "light" and "heavy" soliton. The particular type of the soliton 
is governed by dimensionless parameter: 
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FIG. 1. (Color online) The phase diagram indicating the different propagation regimes of soliton. The shadowed regions (a) 
and (b) correspond to fast light and fast heavy solitons correspondingly, while the dashed region (c) pertains to heavy slow 
solitons. The upper hyperbola, marks the boundary of the region of parameters where a soliton can pass over the impurity, see 
Eq.© 



a. 



2a 

v 



(5) 



We name a soliton the amplitude of which is mush less than its velocity, a « 1, a light soliton. Correspondingly heavy 
solitons are those for which the inverse inequality is satisfied, a ^> 1. This difference can be seen in FigJT] the region 
marked as (a) corresponds to the light solitons while regions (b) and (c) correspond to the heavy ones. The properties 
of light solitons are very close to those of a linear wave packet: almost all of its (positive) energy is contained in the 
first term in the r.h.s. of Eq.@ (i.e. in its kinetic energy). On the contrary, a heavy soliton behaves mostly like a 
classical particle and its (negative) energy is concentrated mainly in the second term in the r.h.s. of Eq.Q. It means 
that the energy of the nonlinear interaction of the quasiparticlcs bound in the heavy soliton significantly exceeds their 
kinetic energy, and the opposite is true for the light solitons. 

Consider now a soliton incident at t — > — oo from the right (i.e. from x = +oo) and characterized by the amplitude 
a, velocity v (or equivalently by the number of quasiparticles N s and energy E s Q), and also by the phases <po,i- As 
noted above, in this paper we will deal with the case of weak scatterers where the dimensionless scatterer intensity is 
small, and the ratio s/v <C 1 is the main small parameter of the problem. Then one can resort to perturbation theory 
[t| |27J to describe the change of the soliton parameters after each act of scattering. The exhaustive perturbative study 
of the soliton transmission over a single delta-scatterer was carried out in Rets. [l(J. It was shown that in the case 
of attractive scatterer, e < 0, the soliton becomes trapped in an effective potential well created by such an impurity 
and then experiences oscillatory motion in the vicinity of the scatterer. In the case of sufficiently strong repulsive 
scatterer, when the two inequalities e > and v < 2ea are fulfilled, the soliton is always reflected by the scatterer 
(the region of parameters corresponds to the upper right sector in Fig. [T|). The soliton can pass through the scatterer 
only in the case of comparatively weak repulsive interaction, and in our current study we assume that the following 
inequalities hold: 



< 2ea < v. (6) 

We distinguish two cases: i) a fast soliton, with parameters falling into the two marked regions, (a) and (b) in Fig|TJ 
where the strong version of inequality ([5]) is valid: 

ea <C v, (7) 

and ii) a slow soliton, where together with Eq. we have ea ~ v (region (c) showed schematically by dashed lines 
in Fig[TJ. As one can easily see, in the case of weak scatterers, s -C 1, the light soliton is always fast, while the heavy 
soliton can be either fast (region (b) or slow (region (c)). 
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In the first order with respect to dimcnsionless scatterer intensity £, at t — > +00 the amplitude, velocity, the number 
of particles and energy of the soliton remain unchanged. The only changes are in the phases: tp' 1 = 930,1 + O (e). Here 
and further on the primes will denote soliton parameters after scattering. These changes can be readily calculated 
[Toj and as a result the position of the soliton center is asymptotically shifted back as compared to the unperturbed 
propagation of the soliton with the unchanged constant velocity v. The overall value of the backward coordinate shift 
(in the adiabatic approximation) is given by 



d^2e/v 2 . (8) 

For a slow soliton the additional shift of soliton position d is of the order of its width l n ~ 1/a, and for a fast soliton 
this shift d <C ^o- 

Since in the first order in e (adiabatic approximation) the amplitude a and the velocity v of the soliton do not 
change, its energy E s and the number of bound particles N s do not change either (see Eq.Q) and there is no emission 
of quasiparticles from the soliton. But in the second order in e (or rather e/v as will be seen later), the solution at 
t — > +00 represents the transmitted soliton together with a number of quasilinear excitations carrying the number of 
waves N and energy E. The amplitude and velocity of a scattered soliton do change now: 

a' = a+0 (e 2 /v 2 ) , v' = v + O (e 2 /v 2 ) , 

and so do the number of bound quasiparticles, N' s = N s + O (£ 2 /u 2 ), and the energy, E' s = E s + O (e 2 /v 2 ). All the 
parameters describing the incident and scattered solution are related via the conservation laws for total number of 
quasiparticles and the total energy: 



E °= E ° + E > (9) 
N s =N' S + N. w 

In the case of slow soliton the spectra of emitted quasiparticles can only be obtained numerically. However it 
appears that such a limit (when the position shift is of the order of the soliton width) is not a physical one. On 
the contrary, the case of fast soliton admits a detailed analytical description. Here the problem of emission of linear 
excitations can be solved perturbatively 0, [l3[ . In the second order of the perturbation theory the number of emitted 
quasiparticles and their energy have been calculated explicitly: the corresponding results in the cases of heavy and 
light soliton look rather different 

E3- As was mentioned above, a heavy soliton behaves like a classical particle and 
the total number and energy of quasiparticles emitted by such a soliton are exponentially small: 



A^^Q M/ VW«, £ ,AQ li/a e -«>, a »v. (10) 

The behavior of a light soliton more or less mimics that of a linear wave packet. In the leading approximation in a, 
the expressions for the number and energy of quasiparticles emitted by a light soliton are as follows: 



iV~2a 



E 



2 
e a 



a <C v. 



(11) 



B. Transformations of the amplitude and velocity for the fast soliton. 

By virtue of Eqs. © - ([TTj) one can express the energy change, 5E S , and the change of the number of quasiparticles, 
SN S , via their emitted values. However the natural kinematic characteristics of the fast soliton are its amplitude and 
velocity. Now by means of the results of previous subsection we can obtain the mapping for the fast soliton amplitude 
and velocity during an act of scattering on a single defect. According to Eqs. ^ the small changes of the soliton 
parameters after the scattering are given by (see Ref . fl3| ) : 

v 2 

SN S = -N = 2 5a, 5E S = -E = — Sa + av Sv - 2a 2 da. 
These relations determine the relative change of soliton amplitude and velocity 

Sa N Sv E N aN 

— = -*-' "7 = 2+ A 2~> ( 12 ) 

a la v av^ 4a v z 
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yielding the following evident velocity and amplitude transformation rules: 



1 - 



E 



N 
4^ 



1 - 



4a' 



(13) 



N 
2« 



(14) 



see Ref . fl4j . A continuous version of these equations was obtained in Refs.fl3l. Il6j. 

In the case of fast heavy soliton, where the inequalities are ea <C v 2 and a ^> v, we can use the explicit expressions 
for the number N and the energy E of emitted quasiparticles (fTU|) mentioned in the previous Section III Al Here the 
last term in the second of Eqs. (|12p dominates and the transformation relations become: 



1 -2tt 



u » ' -ira/v 



2 , n \ 7/2 / 

-ira/v 



a > u. 



(15) 



The case of light soliton, a -C w, is more subtle. In the leading approximation in a/w from Eqs. (jlip . due to the 
cancelation of contributions from energy and from number of emitted quasiparticles, the velocity does not change at 
all. Therefore the next terms in expansions of Eqs. (|11[) should be taken into account. Corresponding calculations 
are performed in Appendix \K\ and the eventual results are the following: 



N ~ 2a 



1+0 



E 



2 

e a 



4a 2 „l ar 



a <C v. 



(16) 



For velocity and amplitude transformation in the case of light soliton one obtains: 



v' = v 1 



4eV 
3~ 



a' = a 1 



a <C v. 



(17) 



Generally, accounting for the results above the mapping rule can be represented in a unified way: 



a' = a [1 — G(e, a, v)} , ?/ = v [1 — F(e, a, v)] 



(18) 



where 



F(e,a,t;) = 7-3 G(e,a,v), 



(19) 



and 



-, a < 7;, 



a > u, 



G(e, a, w) 



a <C v, 



(20) 



One can see that in the case of light fast soliton, a v, the relative change of the amplitude is determined only 
by a small parameter (e/v) 2 <C 1, while in the case of a heavy fast soliton, a 3> v, this change contains additional 
exponentially small factor. Also the velocity change in both cases contains an additional parameter (a/v) 2 . Thus the 
velocity change essentially exceeds the amplitude change for a heavy soliton, and is much smaller than the amplitude 
change for a light soliton. Let us note that the results very similar to Eqs. (fT8|) - ([20| but for the different model of 
lengthy segment of random media were obtained in Ref. [15J . 



7 



III. THE MODEL SYSTEM OF RANDOM SCATTERERS 



Let us now introduce disorder into the system and proceed from a single scatterer (defect) to a system of many 
scatterers with random positions and intensities. The NLSE for the system with many scatterers takes the form: 



iu t + u xx + 2\u\ 2 u = u Ek S(x — Xk), k = 1,2, 3, — oo < x < oo. (21) 
fc 

The r.h.s. of Eq. ([21j) describes the influence of the point scatterers with random intensities placed at random 
positions x k . We will consider the statistical properties of these quantities separately, starting with the intensities. 



A. The random intensities 



Let us consider the intensities of the defects as the mutually independent random variables with the common 
probability density function, /5(e), and hence the same two first moments: e and e 2 = e 2 , (the bar denotes averaging 
with the probability density p(s)). We assume here that the distribution of e is not extremely exotic so that e and 
£o are of the same order of magnitude. The squared intensity of the k-th scatterer can be also characterized by a 
dimensionless parameter: 



5 k = % - 1. 

Evidently, the first moment now S k = 0, and all second moments are equal to the same value 5f. = 5$. The parameter 
<5q is of the order of (A/eq) 2 where A is the width of the probability density function. 



B. Spatial distribution of scatterers 



As for the spatial distribution of the defects, we will assume that they are distributed independently and uniformly 
within the segment [0,L] with the mean distance / between the adjacent scatterers. The number n of defects on 
the segment [0,L] is random, and the probability p n to find exactly n defects within the segment, is taken to be 
Poissonian: 



A" 



-A 



n=0 



(22) 



where A = L/l is the average number of defects on the segment. Let us consider the probability density to find exactly 
n scatterers at the points x k , where we ordered the positions: 1 < k < n, < x n < x n -\... < xi < x\ < L. It is 
convenient to introduce the new dimensionless variables: 

zi = L~ X (L - xi), z k = L _1 (x fc _i - x k ), 2<k<n, 
see the scheme given in Fig. 2. Now the aforementioned probability reads as: 



Pn({z} n ) = n\6 1 - z k 
\ fc=i / 



n °? n 

Y[0(z k ), Pn({z} n )Y[dz k = l 

k=l n fe=l 



{z} n = (zi,z 2 , ...,z n ). 



C. The averages 



The dynamical quantities of the problem are described by various functions which depend on the number n, set of 
positions {z}, and intensities {S} of the scatterers. In this subsection we will present a simple and universal method 
of calculating the arbitrary averages over the intensities and Poisson-distributed positions of the scatterers. Again, 
the averaging with respect to intensities of the scatterers is denoted by an overlinc: 



8 



X n X„-, X„- 2 X 3 X 2 X, '- 

I 9 O — — — — 9 9 9 1 



FIG. 2. The scheme of a classical particle scattering on a disordered segment. 



/(n, {z}, {6}) = J Pn ({5})f(n, {z}, {5})D n {5}, p n ({8}) = p($i) ■ ■ ■ p(S n ), D n {6} = d6 1 ... dd n . 

Now, to distinguish between the intensity and position averaging the symbol (. . .)„ will be used for canonical 
averaging with respect to positions of the scatterers: 



(f(n, {z}, {8})) n = J Pn ({z})f(n, {z}, {5})D n {z}, D n {z} = d Zl ... dz n . 

By "canonical" averaging we mean the averaging over all the realizations where exactly n scatterers occur inside the 
segment [0,1/]. Using the analogy with the notions of statistical mechanics we can introduce the "grand canonical 
averaging" over all the realizations having different number of scatterers with the weights given by Eq. (|22|) . The 
grand canonical averaging will be denoted by a single pair of angular brackets without the subscript n: 

(f(n,{z},{S})) = J2p n (f(n,{z},{S})) n . 

n 

Double angular brackets will denote the total averaging: 

«/(n, {z}, {5}))) = </(n, {*},{«*})> = {*},{*})>»■ (23) 

n 

In what follows it is useful to recall the expressions for the first two moments of the Poissonian distribution: 

oo oo 

(n) =^2np n = A; (n 2 ) = n 2 p n = A 2 + A. 

n=0 n=0 

We can also obtain the first two canonical moments of the distance between the scatterers: 



( z k)n = J z lPn ({z})D n {z} = — j-p (zl) n = J z 2 p n ({z})D n {z} = j- - 1 2 - ^ , 

(ZjZ k ) n = I z 1 z 2 p n ({z})D n {z} = - — — - — — , j ^ k. (24) 
J (n + l)(n + 2) 

To obtain the latter quantities we have used the following trick. One starts with the well known integral representation 
for the theta-function, 

8(u) = -L f <^t e ^u 
K ' 2ni J c k 

(the contour C in the complex plane k is the line Im/t = —0). Then it is possible to interchange the order of the 
contour integration and the integration over all Zk ■ Note that this trick comes handy for calculating arbitrary canonical 
averages, not only the first moments. 
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IV. THE TRANSMISSION THROUGH A DISORDERED SEGMENT 



A. Dynamics 

In this section we consider the transmission time of a soliton passing through a disordered segment and calculate 
the mean value and variance of the transmission time shift. It is important to mention that in the following analysis 
we account for the linear radiation only when applying the transformation formulas, Eqs. (|18[) . However we neglect 
all the secondary effects relevant to the multiply reflections of the radiation emitted previously by the soliton and the 
weak recurrent action of these waives on the overall soliton dynamics. 

Let Vi > be the initial velocity of the soliton incident from the right on the segment [0, L}. In the absence of 
scatterers the transmission time would be T = L/v\. The inclusion of random scatterers brings about two effects. 
First, the effective distance between the adjacent scatterers increases because of the backward shift of the soliton 
position. We denote such a backward position shift at the fc-th scatterer as d^. Second, a small deceleration of the 
soliton occurs after each act of scattering, which also increases the transmission time. Let Vk > 0, k < n be the 
velocity of soliton incident (from the right) on the fc-th scatterer, and v n +i - the velocity of the soliton after passing 
through the last n-th scatterer (note that the sequence of velocities {vk} monotonically decreases, Vk+i < ffc). Since 
the soliton propagates from right to left (i.e. from +oo to — oo), we assume that Vk denotes the absolute value of the 
velocity. 

The total soliton transmission time T n through the disordered segment containing n scatterers is: 



L — xi , x 1 -x 2 + d 1 x n -i—x n + d n -i , x 



T Q + AT nl + AT„ 2 , (25) 

Vl V 2 V n V n+1 

where the meaning of the quantities AT„i and AT„2 is following. By AT„i we have designated the quantity: 

n , 

AT nl =V^-. (26) 

So it is the contribution to the transmission time delay related to the backward position shift. This delay, AT„i, 
according to Eq. ([8} can be expressed as: 



9c 

AT nl = J2~fi^- (27) 

It follows from Eq. (|I5p and Eq. (|17|) that (ufc+i — ffc)/wfe = (efc/ffc) /(a^), so that in the leading approximation in e 
one can substitute v\ for all the parameters Vk- Then we arrive at the expression: 

2 " 

AT nl ~— Ve fc . (28) 

and the final result for the first order contribution to the delay time is: 

A7^=^. (29) 
The last term in the r.h.s. of Eq. (j2"5]) . AT„2, is as follows: 



AT„ 



fe=i 



(30) 



so it describes the delay which occurs due to the soliton deceleration after each act of scattering. It is written in the 
most general form and as such describes the propagation time of an arbitrary classical particle traversing through 
a disordered segment and obeying some given velocity transformation rules. Therefore the following analysis of this 
contribution in the current and following sections is applicable to any classical particle moving according to a given 
monotonically decreasing velocity mapping. Note that the velocity deceleration need not be small for the validity of 
the method - only the positiveness of the velocity is important. We note that the intermediate calculations and final 
results also look simpler in the general form rather than in the form of expansion in the limit of small e/v. The latter 
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limit is considered only in the subsection IIV CI where we obtain and analyze the mean value and the variance of the 
transmission time delay for a fast soliton. 

In the dimensionless variables Eq. ([30| can be recast as: 




The statistical analysis of the shift AT„2 is more complicated than that of AT„i and the rest of the paper is dedicated 
to the former. As we will see below this contribution dominates for long enough segments with the large average 
number of scatterers. Therefore in what follows we will omit the subscript "2" and write simply AT„ implying the 
quantity AT„2 unless specified otherwise. 

The corresponding delay depends on the number n of the scatterers and on their realization, i.e. on the two sets of 
parameters {z} = {z±, . . . , z n } and {6} = {Si, . . . S n } (the latter enters through the velocities i>fc). Indeed according 
to Eqs. (fl~8)l , the change of velocities is described by the set of recurrent relations: 

Vk = v k -i [l - F(ek-i,a k -i,v k ~i)), a k = a k -i [l - G(s k -i, Ojfe_i, v k -i)] , (32) 

which actually gives Vk as a function of an input velocity v\, input amplitude a\ and fluctuations Si of intensities e% 
of all scatterers I = 1,2, k — 1 preceding (from the right) the scatterer with the number k, 

Vk = $(vi,a 1 ;Si,S2,...,5 k -i). (33) 

We emphasize again that in the case considered, the soliton behaves as a classical particle that moves with a 
constant velocity between the scatterers. The corresponding dimensionless transmission time shift, 



T n({z},{S}) = — ^ = //„,i - VfriZ t , (34) 

k=l 

is expressed via the natural dynamic variables: 

Mn.fc = — — > 0, fc = l,2, ...,n, (35) 

V n +1 Vk 

which arc nothing else but dimensionless shifts of the inverse velocities after passing through the last n — k + 1 
scatterers. Thus we see that the fluctuations of intensities and positions of the scatterers in Eq. (|34p are decoupled: 
accounting for Eq. Q33[) the former enter only through the variables fi n ^k while the geometric disorder enters directly 
through the dimensionless distances Zk between the adjacent scatterers. 

Because the velocities form a monotonically decreasing sequence, the sequence ii n ,k also decreases: [i n ,k+\ < Mn.fe- 
From Eq. (|3"4")l it follows that the shift r„ does not exceed the value /J, n ,i- 

< t 71 ({z},{5}) < (i n>1 . 

The maximum of the transmission time occurs in the configuration where all the scatterers are concentrated at the 
point of incidence (all Xk = L) so that the soliton always moves with the minimal velocity v n+ \. The minimal shift 
equals zero and corresponds to the opposite configuration where all the scatterers are concentrated at the farther end 
of the segment (all Xk = 0) and the soliton always moves with its initial (maximal) velocity V\. 

B. The mean value and the variance of the transmission time shift 

The statistical properties of the dimensionless transmission time shift, Eq. (|34[) . are the main subject of our paper. 
We will show that within the framework of the classical "particle" model one can construct a complete statistical 
description of the transmission time delay and obtain a general expression for its probability density function. This 
problem will be considered in the next Section IVl 

However, sometimes for practical applications one may want to know only the first two moments of the delay. In 
this subsection we will calculate separately the mean shift of the transmission time and its variance. To calculate the 
variance, besides Eq. (|34[) we will also need the expression for the square of the shift in terms of [i n ,k and Zk- 
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Tn (-{>}> {$}) = Mn,l " 2 M«,1 + X! Mn,fe Z fc + E Mn.fcMn.m^k* 

k—1 k—1 k^m 

The canonical averaging of Eqs. f[5^|) , 

^ = (r„({z}, {«}))„, r( 2 ) = («})>„, 
by virtue of Eqs. (j2~4")l leads to the following general expressions: 



(36) 



t^ 5 = Hn,l 



1 ™ 

— y 



k=l 



r (2) - u 2 

'n — H>nA 



k=l 



(n + l)(n + 2) 



fe=l 



Ofc=l 



Eventually the total averages can be obtained by applying the averaging formula, Eq. (|23[) : 



rW = ((r„(M, {*}))) = E^ r n X) 



Tl = 



r (2) = ((^(M,W))) = £p„r^ 



2) 



n=0 



C. Weak scatterers limit 
Eqs. (f3"T)) . (|3"8")l were obtained for the general type of velocity mapping: 



(37) 
(38) 



(39) 



Uj+i = u j C 1 - ^Oj)) 
Here we consider the weak scattering limit where the parameter 



(40) 



C = = F{vi) 

Vl 



(41) 



is small. More precisely, to avoid the error accumulation it should be much less than inverse number of scatterers on 
the segment 



(42) 



This equation serves as definition of the weak scattering limit. Turning to the specific problem of the soliton scattering, 
the function F(v\) coincides with the function F(eo,ax,vi) from Eq. (fT9|) . 

In the first order in £ the mapping relations defined by Eqs. (|32[) . (|35[) become much simpler. For the fast soliton 
they read: 



fc-i 



l — Clfe-i+E^" 

V rn=l 
/ n N 

fj-n.k = C [n- k + l+ ^2 S m 



(43) 



Note that the velocity mapping decouples from the amplitude mapping in this approximation. The canonical averages 
can now be straightforwardly calculated: 



r«=C- t& =?(-)' 

' n S cy i ' n \ 2 / 



(44) 
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One notes that in the case of large number of scatterers «,(«(! 1 < 1, both dimcnsionlcss shift of the transmission 
time, t„ , and its relative standard deviation, 



.(2) 



\0 



(1) 



3/v 



(45) 



are small. The smallness of fluctuations is provided by large number of scatterers n only, while the small deviations 
from the unperturbed transmission time require both the weak strengths of the scatterers e and/or the large initial 
velocity v\ (recall that £q/v\ <C C 1). 

After the next averaging over the number of scatterers, Eq. (|39j) . we obtain: 



oo 



- (1) = E^ 1) = -- 

n=0 



2 



n=0 



3 v 



(46) 



These results are qualitatively the same as those for the canonical ensemble, cf. Eqs. (|44|) . The difference is that here 
the average number of scatterers, A = L/l, stands for n from Eqs.(|44j). and some of the numerical coefficients have 
changed. Finally, going back to the dimensional variables we obtain the weak scattering expressions for the moan 

transmission time shift (AT) , and its standard deviation 5T = ( ^(AT) 2 ^ — (AT) 2 ) 1 ^ 2 : 



(AT)„ = ^T 0) 



ST n = (AT), 



l + 45 2 



3n 



(47) 



for the canonical ensemble, and 



((AT)) = ^T , 6T=((AT))yJ Al{1 + L S2o) , (48) 

for t he g rand canonical ensemble. The expressions above present an improved version of those given in our earlier 
Ref. [26[ where the definition of parameter £ for the case of a light soliton was erroneous. 

Note that the values of small parameter ( for light and heavy solitons are different. Therefore within the canonical 
ensemble, the mean shift of the transmission time is: 



2 nL el 



light soliton a -C v, 

3 v v \vJ 

<AT)„ = ^ (49) 



nL £n /a\ n / 2 / a\ . .. 

ir ? I — exp — 7T— , heavy sohton a>t. 

v v V v ' V v J 



Here and in the two following Eqs. (|50j) and f|5 1[) . we set for simplicity a = a\ and v = v\ as the designations for the 
input amplitude and velocity values. The results for grand canonical ensemble are obtained by replacing the number 
of scatterers n by its mean value L/l. The standard deviation is then found from Eqs. (|T7|) . (@5|. 

Recall now that the results above are relevant for the second order contribution in e, AT„2 (which is due to the 
velocity shifts and deceleration of the soliton) to the transmission time delay which coexists with the first order 
contribution T„i, Eq. (|29[) . (which arises due to the position shifts). Each of the two is dominant in its own interval 
of the segment lengthes. The crossover length, L c , is defined by the expression: 

1 /v\ 2 



light soliton a <C v, 

so \aJ 

L c ~ { (50) 



1 /U\"/2 



exp ^7r— J , heavy soliton a 3> v, 



For short segments, L -C L c , the transmission time delay is mostly determined by the first order contribution due to 
backward shift of the soliton center after each scattering, see Eq. (f29]) . On the contrary, for the long segments, with 
L ^> L c , the velocity deceleration plays crucial role and the delay is described by Eq. (|49|) . Finally, combining the 
results in the two aforementioned regimes we can write the mean transmission time delay for canonical ensemble as 
a general expression: 
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max 



(AT)„ = { 



2ns 



3v 3 



max 



— , — l-j CX P I-7T- 



11/2 



light soliton a <C v, 
heavy soliton a ^> v. 



(51) 



V. STATISTICAL PROPERTIES OF THE TRANSMISSION TIME DELAY 



A. General formulas 



In the previous section the two first moments of the time shift were calculated both in general situation and in 
the weak scatterer approximation. The exact results were applicable to any classical "particle" evolving according 
the prescribed velocity mapping between the scatterers. Here we perform the analysis a step further and obtain the 
probability density function for the soliton transmission shift of a scattered "particle" . 

The latter can be written as a grand canonical averaging of the corresponding delta function: 

p(r) = {{S(r-r n [{z},{5}}))). (52) 
Expressed in terms of canonical averages this probability density is: 



P( T ) = ^PnPnir). (53) 

Here po(r) = 6(t), and for all n > 1 



p n (r) = p{T;n,{5}), (54) 

with the partial probability density 

p(r;n,{6}) = (6(T-T n ({z},{6}))) n . 
The probability density will be explicitly calculated in the following subsections. 



B. Probability density function 



The intensities of the scatterers enter the partial probability density, p(r; n, {5}), through the set of parameters 
fXn t k- This dependence should be taken into account explicitly only during the averaging over all realizations of the 
set {5}. In what follows we will consider the canonical configurational averaging (i.e averaging over the positions of 
the scatterers, z^) and omit the symbol {5} in the argument of the partial probability density. 

For n = 1 after straightforward integration over Z\ one gets: 



p(r;l) = 0(^i,i -t). (55) 

For all n > 2 we also start with the integration over z\. Due to the presence of delta function, the only point 
contributing to the integral is 



t 1 " 

Zl = Zl({z}) =1 ^ z kPn.k, 

which leads to the result: 



Mn,l Pn,l k=2 
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p(r;n) 



| pOC poo I n \ n 

— / ... / 6{zr{z}) e[l-z 1 {z}-y j z k \T\dz k . 



(56) 



Using the integral representation for the 6— function from subsection IIII C[ we can rewrite the partial probability 
density as: 



p(r;n) 



i\ 1 



dux f d,K2 



A*r*,l (27Tl) 2 J Cl Ki Jc 2 K 2 

n 

n 



fc=2 







| J dz k exp 









exp 



JKl ( 1 — 



+ «2 1 



All integrals over {z} converge because the sequence \x n i decreases. After integration we obtain: 



p(t;ti) 



i\ 1 



dni f dtt2 



fi n ,l (2m) 2 J Cl Ki Jc 2 K 2 



■ exp 



IKl ( 1 J + IK2 



n 

fc=2 



Pn,k , . / -, P"n,k 
JKl 1 h IK2 I 1 — 



Mn,l 



The next step is the integration over k,\ and K2- We start with integration over k\. The integrand has simple poles 
at points: 



K\ = K2 1 - 



A*ri,l 



, k — 1, 2, . . . , n, 



which lie at the upper half plane of k±. If r > p, Ut i, we close contour C\ in the lower half-plane and the integral is 
equal to zero. In the opposite case, r < p, n i, the integral is proportional to the sum of residues in all the poles: 



n) = — i 2^ 



n 



dfi2 { IK2 , , 

- — — exp<^ r - (//„i - OU > n >2- 



fe=l 



The integral above differs from zero only for r > (/i n ,i — fJ> n ,j) ■ in this case one should calculate the contribution from 
a single pole K2 = of the order n. The result is: 



n>2, 



where 



/j(t; n) = 0(/i„,i - t)0(t - (/i„,i - /i TlJ )) 



n[r - (/i„j - p. n j)] r 



fJ"n,j Jo 

The specific structure of Eq. (|57p enables us to represent it as a ratio of two determinants: 



fj(i",n)dT = ffi/- 



(57) 



(58) 



p(r;n) 



W n {p-n} ' 

where W n {p, n } is the Vandermonde determinant based on the powers of p. n ,k- 



(59) 



W n {p, n } 



u n ~ 2 


Li"' 1 
Pn,2 ■ 

n-2 

Pn,2 


.n-1 

• ■ fn.n 

• ■ r"n,n 


P-n.l 


M",2 ■ 


• • P>n,n 



1 1 



(60) 
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and W 7 ({fi n } is obtained from the expression above by replacing the first line with the corresponding functions fk 
given by Eq. (J55J : 



fl J2 

n-2 n-2 



1 1 



Indeed, ordering the differences in the denominators in the l.h.s 
summing up all terms, we get the common denominator 



fn 
H>n,n 



1 



(61) 



of Eq. ([59l) . so that all of them are positive, and 



ipn,k ~ Mn,j) = W n {fJ,„}, 

l<k<j<n 

while the corresponding numerator is nothing but the expansion of W£{fJ, n } with respect to its first line. Thus the 
delay probability density is defined by Eqs. (|53l) , (|54[) with partial probability density function p(r; n, {S}) given by 
Eq.©. 

We recall that the scaled transmission time delay, r, enters into this formula through the function /, see Eq. (|58l) . 
while the realization of scatterer intensities, {<5}, enters through the variables p which are determined via soliton 
velocities Vk obeying the velocity transformation, Eq. (|33|) . Eqs. (|22|) and (|58|) automatically provide normalization 
condition of the probability density of the transmission time shift: 



p{t)cLt = 1. 



In particular case of identical weak scatterers, < £j = e <C 1, the results obtained look essentially simpler. Indeed, 
in this case according to Eq. (|43| we get: 



fJ>n,k = C( n -k + 1), 

and in terms of the scaled time shift, f = t/£, the partial probability density p(f,n) = C,p{r, n) reads: 



(62) 



p{r\n) = < 



r 

l-r), 



n = 0, 
n = 1, 



(63) 



0(n-f)^f-j + l)P(f;n,j), « > 2- 



Here all P(f; n, j) are (n — 1)— power polynomial functions of f : 



P(r;n,j) 



n (f-j + iy 



n - j + 1 



(64) 



n ( fc -^) 



All partial PDFs, Eqs. ((63j), arc universal. The fact that we consider either light or heavy soliton, as well as the 
dependence on a particular set of initial parameters a\, v\ and £o, is reflected only in a particular value of scaling 
parameter £, Eq. (j4T|) . Therefore in the case of identical weak scatterers the canonical PDFs given by Eqs. (|63|) . not 
only describe the statistical properties of the time shift of an NLSE soliton but are also applicable to any classical 
object with the velocity transformation law 



v j+1 = v j (l-0- (65) 

With the growth of the partition number n, the singularities in PDF become weaker. The zeroth term, n = 0, 
contains a delta-singularity at the origin, the first term has a jump at f = 1. All higher partial PDFs, p(r;n), with 
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numbers n = 2, 3, ... arc continuous together with their first n — 2 derivatives. However already n — 1— th derivatives 
have a jump at the points f = 1,2, n. For the points f = 1, 2, n — 1 this statement is evident (see Eqs. (1631) . 
([64]) '). For the point f = n it immediately follows from the representation given by Eqs. ([59)l . ([60]) . 

In the general case, performing the summation of canonical PDFs. Eq.(|59j). with the Poissonian weights given 
by Eq. (|22p . we arrive at the grand canonical PDF, p(r), for the transmission time. In the same approximation 
of identical weak scatterers the grand canonical PDF, p(f), (in scaled variable f) depends on a sole dimensionless 
parameter A = L/l which is nothing but the mean number of scatterers in the segment. In Fig|3]wc displayed the 
truncated grand canonical PDF, 

oo 

P(r) = Y,PnP( f > n ) = P^) - e ~ H ( f )> ( 66 ) 

for A = 10, i.e. we extracted the singular contribution at the origin. The truncated (regular) PDF is continuous on 
a whole semiaxis f > except for the point f = 1. At this point the PDF has a jump with the magnitude which 
according to Eq. (|2"2l is equal to 0.000454. This small jump is shown in the inset panel of Fig. [31 

The results above allows one to gather all statistical information about the time delay in the disordered segment. 
One can verify directly (see Appendix [B]) that the expression for PDF given by Eq. (f5T)|) leads to the same expressions 
for the two first moments of time delay as was given by Eq. (f3"Tj) before. 

C. The cumulative distribution function 

Another important quantity of interest is the cumulative distribution function, i.e. the probability P(tq) that the 
dimensionless transmission time shift, r = AT/Tq = t>iAT/L, does not exceed a fixed value tq. This probability, 

P(t ) = P(t<to) = ((6(t -t))), 

is equal to 

oo 

P(t ) =PQ + ^PrK, (67) 

71=1 

where 

0„ = K{tq,{6}) = (tf„(r - T n [{z},{S}])) n . (68) 
Direct calculations similar to those used in the previous subsection (see Appendix [C| lead to the result: 



P(f) 




FIG. 3. Truncated grand canonical PDF for the scaled transmission time delay. The jump of the PDF at f = 1 (magnified) is 
shown in the inset. 
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fc=l 

where 



5j( r o) = f(r - 0„,i - /u n ,j)) . (70) 

fJ"n,j 

In Eq. (|69| we used the same notation for W% as in Eqs. (|57|) - ([59)) . The scaled transmission time delay f enters into 
Eg.(|M]) both explicitly and through the functions g, Eq. tjTO")) . while the realization of scatterer intensities {5} enters 
by means of variables /i, which, in turn, are determined via the soliton velocities Vk obeying the velocity mapping, 
Eq. (J33l). 

Another way of obtaining this probability is the direct integration of the probability density function: 



tfn(To) = / Pn(T)dT. 

Jo 

The calculation of this integral is reduced to the (straightforward) integration of the function fj . The latter is equal 
to: 



fj(T;n)dT = 6(t q - p n ,i)Pn/ + ^(Mn,i - 7o)Sj(to), 



and results in the expression: 



V n {T ) = 0{T - p n .i) + &(/in,l ~ r o)-fTT-f 7" = ^(^0 - AVi) 1 - 7" + TTT 7 r, (71) 

which at first sight looks different from Eq. (|69|) obtained above. However simple calculations (see Appendix [D]) 
confirm the equivalence of these two formulas. 



VI. CONCLUSION 



In this paper we have studied the propagation of the envelope NLSE soliton through a segment containing weak 
point-like scatterers. Both the positions of scatterers and their intensities were assumed to be random and the 
concentration of the scatterers was assumed to be small (the mean distance between the scatterers is much larger 
than the soliton width). 

For a relatively short segment, the transmission time delay of a fast soliton is mostly determined by the shifts of the 
soliton center after each act of scattering. However for sufficiently long segments the main contribution to the delay 
stems from the shifts of the amplitude and velocity after each scatterer. The crossover lengths separating relatively 
short segment from a sufficiently long one have also been obtained for both cases of light and heavy solitons. 

We have developed a method for calculating the statistical properties of the transmission delay time. This method 
is applicable not only to the particular problem of the NLSE soliton transmission but also to the problem of forward 
scattering of an arbitrary classical particle provided that the velocity change during an individual act of scattering is 
known analytically. 

The exact probability density function of the soliton transmission time delay and its two first moments have been 
found. In the case of identical scatterers we have obtained the PDF as a universal function that depends on a sole 
parameter - mean number of scatterers in a segment. Thus the first term of the series contains a delta-function 
singularity, the second has a jump at some point, and each next term contains a jump in the corresponding higher 
derivative (first, second etc.). The relative fluctuations of the delay time become negligibly small when the segment 
length grows. 



18 



ACKNOWLEDGMENTS 



We are thankful to G.R. Belitski for helpful discussions and M.M. Bogdan for the important criticism. This work was 
partially supported by a joint scientific project No.24-02-a Of NAS of Ukraine and RFBR, Israel Science Foundation 
(Grant # 944/05), the joint French-Ukrainian project in the framework of scientific cooperation between NASU and 
the CNRS. SD, JP and AK would also like to acknowledge the support from the UK Royal Society. SG would like to 
acknowledge the support obtained from School of Engineering and Applied Science at Aston University, Birmingham, 
UK for his visit to UK. 



Appendix A: The number and energy of waves emitted by a light soliton 

In the second order of the perturbation theory with respect to e, the number N of emitted waves and their energy 
E are proportional to the following integrals (see Refs. [JOS an< ^ a ^ so 03 )■ 



N 



7r e 



F(y,a)dy, 



E 



28 c 



e v y F(y,a)dy, a 



2a 

V 



where F(y,a) is the renormalized spectral density of emitted waves: 



F(y,a) 



cosh 



7r [y 2 - 1 + a 2 
4a 



Let us sketch a scheme of a general asymptotic analysis of the moment integrals of type (|A1[) , 



(Al) 



(A2) 



y n F(y,a)dy, 

i 

for a light soliton a -C 1. In this case the main (power law) contribution comes from the two peaks of the spectral 
density (|A2j) . y± = ±1 + 0(a 2 ), but the contribution from the right peak, y+ ~ 1, is generally smaller by a factor of 
a 4 . In the vicinity of each peak, y±, one can make a substitution: 

x = ^(y 2 +a 2 -l), 

which is of course denned together with its inverse y±(x). After the substitution one can write: 



/oo 
f±(x, a)sech 2 (x/a) dx, 
-DO 

where 



f±(x,a) = y&x) [(y ± (x) - l) 2 + a 2 ] 2 \y' ± (x)\. 

Because the a-dependence of the functions f±{x, a) is weak (algebraic) while the hyperbolic function is exponentially 
localized [a <C 1), the asymptotic expansion of the above integral can be performed using an analogue of Watson 
lemma, i.e. by means of substituting the formal Taylor expansion of f±{x, a) at x = and consequently integrating 
term by term. The result is then becomes as follows: 



^ (2*)! L a ^ (2*)! C2ka ' 

where Co = 2 and the coefficients C2k, k = 1,2,..., arc expressed via Bernoulli numbers B^k'- 



n2k-l _ 2 
C2k = 2 2fc-2 ^ l B2fc l ■ 

Substituting the corresponding even derivatives of functions f± evaluated at x = and developing them in series in 
powers of a we obtain the sought asymptotic expansions for each moment n of the spectrum. The results for n = 0, 2 
yield formulas (|TB|) . 
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Appendix B: The calculation of moments 

Consider first a fixed configuration of the scatterer intensities. The functions fj(r), Eq. (|58|) . introduced in subsection 
IV Bl divided by /x™" 1 become the probability densities themselves. The two first moments of r calculated with these 
probability densities are: 



ra-l , 



1 



2 ^(-^n,lf^n,j 



n + l VnjJo n + 1 (n + l)(ra + 2) 



Canonical averages r^ 1 ', and t^ 2 ' 1 , calculated with the help of canonical probability density p(r,n) (see Eq. (|57|) or 
Eq. (|59l)) are: 



(1) _ (2) _ 2 _2^iW^{M 2W^ j {M 



r(2)j 



n + 1 ' 



n+1 (n + l)(n + 2) : 



(Bl) 



where Wn {//„} is obtained from W„{[i n } = Wn {/"n} 5 Eq. (|60|) . by replacing all /i™ j, 1 in the upper row with jU™"^ m 1 . 

All the generalized Vandermonde determinants WJ, are proportional to the initial one W n {fi„}, Eq. (|M| . 

Corresponding coefficients are expressed as symmetric polynomials of fj,„ &: 



n 



fe=i 



\k=l / l<j<k<n 



Substituting these expressions into Eqs. (|B1[) after some straightforward algebra one arrives exactly to Eqs. ([57)1 
obtained in the body of the paper. The additional averaging over {6} completes the calculation. 



Appendix C: The cumulative probability distribution 

Explicit form of the partial probability see Eq. (p5|) . is: 

/>OG poo I °° \ / n 

Jo Jo \ k=i J \ k=i 

With the help of the integral representation of ^-function it can be written as: 



.kZk 



■0,. 



r r dK2 e ^ 2(T0 _ Mr(il) -pr r 

Jd «i Vc 2 K 2 J- =1 7o 



s iZfc(^„,fcK 2 -Kl)^ 



To ensure convergence of all the integrals over Zfe, we choose contours C; (i=l,2) so that not only Imft^ < but also 
Im(ni — Hn t kK-2) < for all k < n. Then integration over all Zk leads to the expression 



c 2 «2 



^[(Kl - Hn,kH2) 



Closing contour of integration of the first the internal integrals through the upper half plane n\ and taking into 
account that all the poles, k,\ = 0,n n ,i K 2, ■ ■ ■ > "n,n K 2i- be within the closed contour Ci, we obtain: 



2iri 



\ 



fc^ ^ 



(-iy 



Y[f-ln,k 

fc=l / 
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Here the prime over the first product in brackets means that in the pole (where k = j) the second term in the 
denominator, fj. n k, must be omitted. As results for i9„ we obtain: 



n n\ f d,Ki 



2^ Jo, 4 +1 



k=l 



a i«2 (tq - Mre , 1 +f"n , j ) 



j = l ^™>2 fc^j ^"'^ 



Now if we close the contour of integration in the upper plane K2 and take into account the sole pole of the (n + l)-th 
order at the origin we arrive at: 



q a , ,(Hn,i-T ) n L ^^(r -(//„,i -/i„,j))(r -(//„,i -/i„,j)) n 
V n = V(To - Hn,l) n + ^ n > 



fc=l 



(CI) 



which is equivalent to Eqs. ([TO")) presented in Subsection IV CI 



Appendix D: Equivalence of the results given by Eqs. (|69|) and (|7ip 

To establish the equivalence of these formulas wc firstly note that the last term in the r.h.s. of Eq. (fTTj) coincides 
with that of Eq. (|55|) . Then the second multiplier in the first term of Eq. (J7TJ) can be presented as: 



W„{At«} W n {//„} : 



(Dl) 



where 



, / \ n-1 ( T - - Mn,j)) T 

^■(To) = /Vj 



The last ratio in Eq. (|D1|) is a n-th power polynomial of r that vanishes together with all its first n — 1 derivatives at 
the point tq = // n ,i- Therefore one gets: 



C(7D - A*»,i) n . 



To find the constant C, we differentiate this identity n times and obtain: 



1 



Recall that Wn"' 3 is defined by Eq. ([6Tjl . The numerator in the last ratio can be easily calculated and equals to 



^Mn, ft 

fc=l 



21 



This result together with two previous equations establishes the identity of the two forms of probability "d n (To) given 
by Eqs.flSH]) and ([71]). 
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